Non-equilibrium dynamics in amorphous SisB 3 N7 
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We present numerical investigations of the dynamics on the energy landscape of a realistic model 
of the high-temperature ceramic a-SiaB3N7. Below a critical temperature T c ~ 2000 K the system is 
no longer in equilibrium, and we predict that the material has a glass transition in this temperature 
range at high pressure. Analyzing the two-time energy correlation function shows aging in this 
system, which is linked to the geometrical properties of the energy landscape. 

PACS numbers: 02.50.r, 05.90.+m, 61.41.+e, 61.43.Fs, 64.70.Pf 
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Introduction: Amorphous nitridic ceramics con- 
taining both silicon and boron such as a-SiBN3C or a- 
SiaBaNv are synthesized via the sol-gel processjjl) and are 
one of the most exciting new classes of high-temperature 
materials. Up to now, no crystalline form of these com- 
pounds is known. Under standard conditions a-SiBN3C 
is thermally stable and amorphous up to « 2100 K, pos- 
sesses excellent mechanical and elastic properties, e.g. a 
high bulk modulus of ca. 200 — 300 GPa, and is sta- 
ble against oxidation in CVatmosphere up to 1700 K|jJ. 
These compounds appear to form covalent networks with 
a homogeneous distribution of the cations at least down 
to a scale of 1 nm[^| and pose many fascinating questions 
regarding their structure and physical properties. 

Probing the glass transition of a-SisP^Nr, the basic 
representative of nitridic ceramics, is experimentally dif- 
ficult, since, under standard conditions, decomposition 
takes place at T w 1900 K, i.e. before the ceramic melts. 
This material is currently considered for high temper- 
ature engine applications and its aging properties are 
therefore of clear technological relevance, beside hav- 
ing theoretical interest within the general framework of 
glassy dynamics. 

Determining whether a possibly non-ergodic system 
has "for all practical purposes" reached thermal (quasi- 
)cquilibrium is not straightforward and possibly con- 
stitutes an ill-posed question. Physical properties of 
amorphous systems are known to drift with the time 
t w , or age, elapsed since the quench into the glassy 
phase. For short observation times t b s <C t w , the drift 
is undetectable and a state of quasi-equilibrium is re- 
vealed by the approximate validity of the fluctuation- 
dissipation theorem. Concomitant to the violations of 
the fluctuation-dissipation theorem for t b s > t w , the 
correlation and response functions acquire an additional 
dependence on t w . This breaking of time translational 
invariance has been observed e.g. in the magnetic suscep- 
tibility of spin glasses, both in experiments |j| and model 
simulations |Q, ||, in measurements of C p for a-SeQ, and 
also in simulations of the dynamical structure factor of 



e.g. a-SiC>2 above the glass transition temperature 0. 

To detect ergodicity breaking we use 1) the specific 
heat CV, which we calculate in three different ways, all 
agreeing in equilibrium but markedly differing if ergodic- 
ity is broken, and 2) the two-time energy-energy average 
fiitw, t b s ; T), and the related two-time autocorrelation 
function CE(t w ,t b s ;T). In quasi-equilibrium, the for- 
mer equals one and the latter equals a generalized stan- 
dard equilibrium specific heat kBT 2 Cv(t w , tabs', T). The 
age dependent Cy has been studied experimentally, e.g. 
for charge-density- wave systems but does not appear 
to have been theoretically explored outside of two-level 
systems at very low temperatures ||. 

Since aging is linked to the complexity of the energy 
landscape of the system, we have investigated some as- 
pects of the latter, emphasizing their relation to the non- 
equilibrium dynamics. 

Model and Techniques: The model of a-SisBsNy 
consisted of 162 Si-atoms, 162 B-atoms and 378 N-atoms, 
respectively, in a 19.1 x 19.1 x 19.1 A 3 cubic box. As an 
interaction potential, we employed a two-body potential 
from the literature [jwj based on an- initio energy calcu- 
lations of hypothetical ternary compounds which repro- 
duces experimental data regarding the structure and vi- 
brational properties of the binary compounds SiaN.4 and 
BN and of molecules containing Si— N— B units. 

The starting configurations for our simulations 
were generated by relaxation from high temperature 
melts JO] . The simulations were performed at fixed tem- 
perature and volume, with a Monte-Carlo algorithm us- 
ing the Metropolis acceptance criterion. In each update, 
an atom is randomly selected for an attempted move in 
a random direction, and with an average size chosen to 
achieve an acceptance rate of w 50 %. One Monte-Carlo 
cycle (MCC) corresponds to N atom = 702 such individual 
moves. Note that the kinetic energy (3/2fcT per atom) 
does not appear in MC-simulations, and that all quanti- 
ties studied relate to the configurational energy. 

The temperatures investigated ranged from 25 to 7000 
K. For each temperature up to 3000 K and above 3000 
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K, 9 and 3 runs, respectively, of length t to tai = 2 x 10 5 
MCC were performed. In addition, for selected tem- 
peratures, ensembles of 100 runs of length t total = 10 6 
MCC were studied. The energy as function of time was 
registered every 10 MCC. Along the individual trajec- 
tories for T = 250, ...,7000 K, halting points xh were 
chosen, from which both conjugate gradient minimiza- 



tions (xh 



and a set of 10 stochastic quenches 



(T = 0K MC-runs) followed by conjugate gradient min- 

imizations [xh — ► xq — ► x m in) were performed. In the 
following, tinit ~ 1000 MCC is the initialization time of 
the MC-simulations needed for the system to reach equi- 
librium in the ergodic regime (i.e. at high temperatures), 
while t w > U n it is the waiting time before the observation 
begin. 

Ergodicity : To investigate ergodicity versus temper- 
ature, we studied the specific heat Cy and the two-time 
energy-energy average 

AU + rp\ {E(t w )E(t w + t b s )) ens /m\ / 1 \ 

<p{t w ,t obs ;T) = j— (T), (1) 

\tLi\tw )ry[t w ))ens 

where the subscript " ens" always denotes an average over 
all trajectories. Cy was calculated using three different 
computational prescriptions. First 



C V (T) 



9{(E) te[t ^ ftotal] )ens(T) 



dT 



(2) 



where time averaging extends from t w to the end of the 
simulation t to tai and the temperature derivative is per- 
formed after the averaging. Secondly 



Cy{T) 



{E(T + AT; t^))tg[t m ,t ro +t ob ,];i oi ,,«^ 
s ' ' ' ~~ 2AT 

(E(T - AT; t w )) te [ tw ^+t obs ];t oba ^:t v 
2AT 



(3) 



This emulates a step experiment where the system ages 
at temperature T. The time averages over the obser- 
vation time tabs "C t w are performed at temperatures 
T± AT, where AT « 0.1T. Finally we gauge the energy 
fluctuations in [t w , t w + t (, s ] by calculating 



Cy(T) 



{{E 2 )te[t w ,t a +t obs ] 



l)ens(T) 



fc B T 2 



(4) 

for a range of observation times t b s which straddles t w . 

C v has no i t, s dependence. By way of contrast, when 
increasing t obs past t w the observed dynamics in C v 
changes from quasi-equilibrium to off-equilibrium (cf. in- 
set in fig.|l|). Cy, which mimicks an experiment per- 
formed after some relatively long equilibration time t w , 
likely yields the most "realistic" value for the specific 
heat for all temperatures. As shown in fig. [l]the above 
prescriptions yield, as expected, almost identical results 
in the high-T ergodic dynamical regime, but differ at low 




2000 4000 6000 8000 10000 

T[K] 

FIG. 1: Temperature dependence of the specific heats Cy b ' c . 
The waiting t w and the observation times tobs were 10 5 MCC 
for C v (D)and C v ; and t w > 10 5 , t obs = 5 ■ 10 3 for C v (o). 
Inset (note the different y-scale): C v for t w = 10 4 , t b s = 10 4 
(o),10 5 (•). Note that C v w C h v for t a 



t ba > tw the two quantities differ. 



1.015 



-* 1.010 



1.005 



1.000 




100 



1000 10000 

^obs 



100000 



1e+06 



FIG. 2: Observation and waiting time dependences of the 
two-time energy-energy average <p(t w ,t bs',T) for 1250 K, for 
an ensemble size of 100 runs (raw data). The inset shows the 
two-time autocorrelation function CE(t w ,t b s ;T — 1250i("). 
Since, even for 100 runs (corresponding to ca. one year of 
CPU time on an AMD 1800+ MP processor), the scatter in 
Ce is relatively large, the data in the inset are averaged over 
ten time steps. 



T. This indicates that below T c w 2000 - 3000 K er- 
godicity is broken. Further evidence stems from the ob- 
servation that for T < T c , the motion is subdiffusive, 
while for T > T c standard diffusion is observed: For 
T > T c w 2100 K, the diffusion coefficients for B, Si 
and N, follow a power law D oc (T — T c ) 7 with 7 = 1.7, 
showing structural freezing-in. Similarly, the relaxation 
times associated with the bond survival probabilities of 
B-N and Si-N bonds display a rapid increase below T c . 

Repeating these investigations for a large number of 
volumes, we find a similar freezing-in of the structure for 
T « T c . Furthermore, we determine a critical point in the 
liquid-gas region of the ternary system (p cr « 0.7 GPa 
and, T cr w 3700K ). Since the tendency to decompose is 
greatly reduced for a supercritical fluid, we predict that 
a-Si3B3N7 should exhibit a glass transition at a temper- 
ature Tg ~ 1700 — 2000 K and at a pressure of pc > 1 
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GPa. Up to now, high-pressure experiments have only 
been performed for T w 1000 K and p w 2.5 GPa. 

For T > T c , the two-point correlation function always 
remains very close to the equilibrium value 1. The ag- 
ing behavior in the glassy phase is shown in fig. ^| for 
T = 1250 K, and for three different waiting times t w = 
3 • 10 3 , 10 4 , 10 5 . In the non-equilibrium regime t a b s > t w , 
4> is seen to deviate strongly from its equilibrium value 
cj> eq = 1. The closely related autocorrelation function 
CE(tw,tobs',T) = (E(t w ) ■ E(t w + t b s )) ens — (E(t w )) ens ■ 
(E(t w + t b s )) ens also exhibits the expected aging behav- 
ior, i.e. a marked decrease to zero from an almost con- 
stant value (oc Cy(t w )) once t Q b s exceeds t w . This mono- 
tonic dependence on t w of the time range t a b s <E [0, t w ] 
during which (quasi-)equilibrium behavior is still ob- 
served, correlates with the stiffening of the response of 
the system characteristic for aging processes: The longer 
the system is allowed to equilibrate, the longer is the sub- 
sequent time range during which equilibrium-like behav- 
ior is observed. This effect concurs with our observation 
that for T < T c we can fit E(t; T) (and also (E(t; T)) ens ) 
over the interval [ti n it,ttotai] as a logarithmically de- 
creasing function, E(t;T) = E {T) - A(T) In (r^rjj ■ 
Neglecting the fluctuations compared to the drift, one 
has (f>(t w ,t b s ;T) sa ^e^")!!^)^ . Expanding 4> for 
tabs < t w then yields cj)(t w .t b s ;T) w 1 + i E u\i ^f- 
Thus, 4>(t w ,t b s )T) substantially deviates from 1 for 
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> t w , as observed in the simulations. The inset 



shows CE(t w ,t b s ;T) plotted as function of the scaled 
variable t b s /t w . As t b s increases, the data appears to 
collapse on a single curve, indicating that t b s /t w scaling 
can be expected to hold asymptotically. 

Energy landscape: Finally, we would like to link 
the non-equilibrium behavior to the properties of the en- 
ergy landscape of a-Si3BaN7. Figure || shows the average 



energy (E(t; T, a^ n )) e n S of local minima x^' m found by 
applying a conjugate gradient algorithm for logarithmi- 
cally spaced halting points along several trajectories as a 
function of time for different temperatures. 

We note that (E(t;T,x^ in )) ens decreases logarith- 
mically with time for T < T c analogously to 
{E(t;T)) ens A fit of the logarithmic slope yields 

A(T) = 76.29 -T— 134.56 • T 2 , which qualitatively agrees 
with the low temperature expansion of {E(t;T)) ens for 
the so-called LS-tree models Q, suggesting that the 
landscape of a-Si3B3N7 might possess some hierarchical 
aspects in that energy range relevant for T < T c . 

For fixed simulation time, the deepest local minima 
are reached for T = 1750 K, which lies right below 
T c . We find a similar behavior for the average en- 

(2) (2) 

ergy {E(T; x min )) of the local minima x min found after 
quenching plus gradient minimization starting from the 
holding points Xh , shown as a function of temperature in 
fig. ^. We clearly recognize a minimum in this curve at 
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FIG. 3: Time dependence of the average energies 
{E(t;T, x^-„))ens of the minima x^] n for selected temper- 
atures T = 750 K, 1250 K, 1500 K, 1750 K. 
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FIG. 4: Temperature dependence of the average energies 
((E(T;x%l))t obs }ens of the minima (Full squares) . The 
open squares are values of reference energies E re f(T) (see 
text). The inset shows the temperature dependence of the 
excess energy A (see text). 



T w 1750 K, and the largest increase occurs at T w T c . 
Analogous observations are well-known from e.g. global 
optimization studies of complex systems, where it has 
been found that reaching the deepest local minima us- 
ing Monte-Carlo-type search algorithms is achieved by 
spending most of the search time in the temperature in- 
terval slightly below the glass transition temperature p3[ . 
Thus this result serves as another confirmation that er- 
godicity breaking is taking place at T T c . 

The second curve in fig. |^ depicts E re f(T) = 
((E(T,x h )) to J ens - (3/2)A^ Qtom /c B T. Note that if the 
motion around the minima were purely harmonic, 
E re f(T) would be the average energy of these minima. 

However, comparing E ref (T) with {{E{T,x y m ' m )) tobs ) ens , 
we note that an excess energy A(T) = E re t(T) — 

((E(T,x^ in )) to J ens > exists. The growth of A(T) 
is monotonic in T and most pronounced for T w T c . 
Clearly, in this temperature region the landscape below 
the halting points, which is "felt" by the walker, changes 
and a substantial shift to higher " reference" values for the 
vibrational contributions to the potential energy occurs. 
Discussion: The computational analysis of a-Si3B3N7 
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using Cy b ' c and (f>(t w ,t i, s ;T) shows that this amor- 
phous material can be expected to exhibit a glass transi- 
tion with a concurrent break in the ergodicity at about 
T c « 2000 K if pressures high enough to prevent de- 
composition are applied. Regarding its structural dy- 
namics and aging properties for T < T c , a-SiaB3N7 ex- 
hibits a general behavior similar to standard test systems 
(Lennard- Jones, a-SiC^) and CDW systems, insofar as 
we observe a freezing-in of the structure, and a waiting- 
time dependence of the two time correlation function and 
specific heat. 

This aging phenomenon is related to the slow non- 
exponential relaxation dynamics on the energy landscape 
for T < T c , resulting in a logarithmic drift towards lower 
energies. This applies both to the actual trajectories and 
the time-sequence of observed local minima. Indepen- 
dent of this aspect of the dynamics, we find that start- 
ing around T c the average potential energy greatly ex- 
ceeds the value associated with harmonic vibrations at 
T c . 'Thermodynamically', this makes itself felt as a peak 
in the specific heat, which is often associated with the en- 
tropy due to an increased availability of additional amor- 
phous configurations. 

One possible origin for this excess energy is trap- 
ping [fl4| [l~5| due to an approximately exponential 
growth in the effective local density of states gi oc (E) oc 
exp(a(E — E m i n )) , with the growth factor a = 1/T tra p ~ 
1/T C . In such a case, the exponentially growing region 
of the deep pockets of the landscape becomes invisible 
for the random walker for T > T trap , and the top of this 
region serves as the reference for the vibrational energy 
contribution. A(T) would thus refer to the depth of the 
exponentially growing pockets with T tra p < T. Studies 
of amorphous networks ]l6[ and polymers |l7[ on lattices 
have suggested that such trapping might contribute to 
the glass transition in structural glasses. 

Alternatively, the phase space volume of the locally 
ergodic regions around the saddle points might substan- 
tially exceed the one associated with local minima, such 
that these saddles serve as reference, as has been sug- 
gested for small Lennard- Jones systems Jlgfl. Here, A(T) 
would correspond to the average energy difference be- 
tween the relevant saddle points and the nearby local 
minima. 

Landscape studies using the threshold algorithm Jl9[, 
which can be used to estimate both local densities of 
states and minima, and the relative size of minimum and 
transition regions |2fJ could resolve this issue, but are not 
yet computationally feasible. Preliminary explorations of 
the region "below" the halting points by performing 10 
quench runs for each halting point x$ concur with the 
results obtained in[^0| and show that a typical halting 
point is associated with only one rather circumscribed 
basin containing many similar minima, and not with a 
large region of the landscape encompassing very different 
structures. The latter would be expected if the dynamics 



were dominated by high-lying saddle points. 
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